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A computer program is described for calculating Bessel functions J n (z) and / M (z), for 2 complex, 
and n a nonnegative integer. The method used is that of backward recursion, with strict control of error, 
and optimum determination of the point at which to begin the recursion. 
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1 . Method 

Given a complex number z and a positive integer NB, BESLCI calculates either 

/„(*), n = 0, 1, . . .,NB-\ 
or 

J n (z), n = 0, 1, . . .,/V£-l 

using double-precision arithmetic. The method used is described in [1 ] l and is based on algorithms 
of Olver [2] and Miller [3], applied to the difference equation 

y«-i= — y w — SIGN -y n +i, (i) 

where SIGN is + 1 for / 's, — 1 for / 's. 

The program sets MAGZ= [|<z|], the integer part of |z|, Pmagz = 0, Pmagz + i = 1, and then 
successively calculates 

Pn+i = SIGN-(-yp„-pn-i), n=MAGZ+l, MAGZ + 2, . . .. (2) 

The sequence is strictly increasing in magnitude [l, sec. 6]. The program takes A to be the least 
n such that \p n \ exceeds a number TEST defined in section 2 and section 3. It then sets 
y^ v) = 0, yj^j — 1/pa, and recurs backward using (1). To the working accuracy, the computed 
sequence y\^\ y\ N \ ... is the recessive solution of (1) which satisfies the boundary condition 
Tmagz— 1. From this solution, the / 's or J 's are found by normalizing: 

Jn{z) = ^ltl « = 0, 1, . . .,/Vfi-l 

I„(z) = ^lfi n = 0,l, . . .,NB-l 
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where 



M= e u ( yi,' v) + 2 J) ( - i) "jf ) /'s, Im 2 >0 

V n=l / 
[A72] 

A* = ></> + 2 J y<y /'s,lmz=0 

/Lt=c" te ( y o v) + 2 2 i*yW J 7 's, Im 2 < 

\ n=l / 

jh= e- ( y<*) + 2 ^ y„ V) ) / '3, Re z > 

/ u=y1 1 v) + 2 l '2 (-l)»y,¥ /'s, Rez = 



M =e- y„*> + 2 ^ (-l)'y„ v) /'s, Rez<0. 

2. Error Bounds for the y n 

For rc > MAGZ, the truncation error in y|, is 

00 STCN Ar ~ r 
r(iv)= v _ v (iv)= n V 

r=v PrPr+l 

see [2]; equations 5.01 and 5.02. This error is bounded by using the following 
Lemma: For n > MAGZ, let 

, _Pn+i_2n Pn-l 

K n - - 



Pn Z p n 

a/id /e£ 



»+i+J^r-i 



X n = r^y 



z > \ z 



(3) 



Le£ p n = mirc (| k n |, X n ). Then for m ^ n, | k m | > p n . 

The lemma, for real z, is lemma 2 of [1]. The proof for complex z is essentially the same. 
The program insures that 



TEST, W2-10 NSIG p# 4 p L+ i, (4) 

where L = max (MAGZ + 1, NB — 1), and NSIG is the maximum number of significant decimal 
digits in a double-precision variable on the computer being used. Then N' is the least n such that 
\p n \> TEST, , and N is the least n ^ N ' such that 



V-T^-TESTV (5) 
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In consequence of (4) and (5), the relative truncation error \T {N) ly n \ is less than \ % 10 -NSI(; for 

all n in the range MAGZ < n ^ /,; see [1, sec. 5]. 

For n ^ MAGZ, it may be impossible to bound the relative truncation error in the above manner, 
owing to loss of precision due to cancellation in (1). Experience indicates that this loss is negligible 
except when the magnitude of v (v) oscillates as n = MAGZ, MAGZ — 1 , . . . , in the back-recursion 

(e.g., calculating 7's, with Re z > Im z). In this case, there will be about D decimals of precision 
in the values > (v) , where D is the number of decimals in ./maczU) (^ma(;zU)) which corresponds 
to NSIG significant figures in the same quantity [1, sec. 5|. 



3. Normalization and Error Bounds for \x 

The equations (3) were chosen to keep cancellation under control. Now 



Jn(z) 



~n r-n 
Jo 



cos (n0)d8 



[4|, 9.1.21. The integrand never exceeds e' Im *' in magnitude, so \J n (z) | ^ e' Im *L Similarly, \/,,(z) \ 
^gl Re z\ m Thus each term of the sums (3) has magnitude less than twice that of the whole sum. 
In fact, these bounds on the magnitudes of J„{z) and /,,{z) are rather weak, and cancellation in 
(3) is less than this would indicate. 

Besides bounding the truncation error of the algorithm, the program provides an estimated 
bound for the truncation error of the normalization sum, defined by eq (3). In the first of these equa- 
tions, this error is 



gov) 



{yo-y^ + 2 J (v„-y<^)). 



For // ^ MAGZ, a bound for the error term y n — y^ is unavailable. For MAGZ < n < /V, 

I Jn ~ y (i K ) | ^ Pnfht I (Pn ~~ 1 ); see [1, sec. 5|. To avoid storing all the />„, the program allows only 
for terms for which n^ N. Here y (A) = 0, and 



Therefore, 



Jn = Pn £ 
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compare section 2 above. Now let 



Then 
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The program sets TEST] 2* 2 • 10 NS1G . The normalization factor p is e "A/magzW [l,sec. 5], so 






./MAGzCz) 



R (N) 



R(N) 



2pl. 



2 P h 1 

(p.V'+l) (p.v-l) 2 |p.v| ~~ (P.V'+I) (P.Y'-D 2 2 



10- 



(6) 



The bound (6) holds for the first, third, fourth and sixth equations of (3); the derivations are 
the same. Similarly, the second and fifth equations yield 






2p 3 v 



(p 2 v-l) 2 2 



i-io- 



These bounds are rather weak, and the error | 5 ( N) lfJL \ turns out to be less than - • 10 _NSI(; . 
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